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The topic of this review is the effects of electron-phonon interaction (EPI) on the transport 
properties of molecular nano-conductors. A nano-conductor connects to two electron leads and two 
phonon leads, possibly at different temperatures or chemical potentials. The EPI appears only 
in the nano-conductor. We focus on its effects on charge and energy transport. We introduce 
three approaches. Eor weak EPI, we use the nonequilibrium Green’s function method to treat it 
perturbatively. We derive the expressions for the charge and heat currents. Eor weak system-lead 
couplings, we use the quantum master equation approach. In both cases, we use a simple single 
level model to study the effects of EPI on the system’s thermoelectric transport properties. It is 
also interesting to look at the effect of currents on the dynamics of the phonon system. Eor this, 
we derive a semi-classical generalized Langevin equation to describe the nano-conductor’s atomic 
dynamics, taking the nonequilibrium electron system, as well as the rest of the atomic degrees of 
freedom as effective baths. We show simple applications of this approach to the problem of energy 
transfer between electrons and phonons. 

PACS numbers: 85.35.Gv, 85.85.-hj, 85.65.-hh, 05.60.Gg, 73.63.-b 


I. INTRODUCTION 

Electron-phonon interaction (EPI) is one of the most 
important many-body interactions in condensed-matter 
and molecular systems^, responsible for a variety of phe¬ 
nomena, from electrical, thermal conduction, supercon¬ 
ductivity to Raman scattering, polaron formation, just to 
list a few^“^. Its effects on the electrical, thermal, and op¬ 
tical properties of bulk semiconductors and metals have 
been intensively studied along with the development of 
many-body theories and experimental techniques. Re¬ 
cent advances in experimental fabrication of meso- and 
nano-scopic structures have generated tremendous efforts 
in understanding the effects of EPI on transport proper¬ 
ties of reduced-dimensional systems^”^^. 

Of special interest are current-induced forces and Joule 
heating in low-dimensional systems, especially in molec¬ 
ular nano-conductors^^“^^. On the one hand, the electri¬ 
cal transport signature of EPI is an invaluable spectro¬ 
scopic tool to study the structural information of molecu¬ 
lar nano-conductors^^’^^. On the other hand, these pro¬ 
cesses are crucial in maintaining the stability of these 
conductors'^, relevant to the continuous scaling down 
of modern electronic devices. Different theoretical ap¬ 
proaches have been developed to study these problems, 
in many cases separately. Recently, it was realized that 
non-conservative nature of current-induced forces pro¬ 
vides an alternative, deterministic way of energy transfer 
between electrons and phonons, or more generally atomic 


motions^^. It is fundamentally different from the stochas¬ 
tic Joule heating. These advances have motivated the 
development of methods treating current-induced forces 
and Joule heating on the same footing^^“^^’^^. 

Equally significantly, there has been an increasing in¬ 
terest in the thermoelectric properties of low dimensional 
systems^^“^^. A starting point of the theoretical treat¬ 
ment is to ignore the effect of EPI, and study the trans¬ 
port of electrons and phonons separately. But how im¬ 
portant the effect of EPI is is a pertinent question, on 
which much of recent work is devoted to^^“^^. Here, we 
will look at this problem using the various approaches we 
have developed. 

EPI is a genuine many-body interaction, the exact 
treatment of which is challenging, if possible at all. One 
natural approach is to perform perturbation calculation 
over a certain small parameter. In the most common 
multi-probe transport setup (see Fig. 1 and Sec. II), this 
small parameter can be chosen according to the strength 
of EPI. This strength can be roughly characterized by the 
ratio between two time scales: the first one corresponds 
to the phonon period, and the second one corresponds 
to the electron dwell time^^ in the nano-conductor. If 
the time electrons spend in the nano-conductor is much 
shorter than the phonon period, the system is in the weak 
EPI regime. The small parameter is the EPI matrix. In 
the other limit, the coupling of the nano-conductor to 
electrodes is the small parameter, over which one can 
perform the perturbation expansion. 
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In this review, we summarize our own effort in devel¬ 
oping and/or utilizing different theoretical approaches to 
study the aforementioned problems in different parame¬ 
ter regimes. We discuss some relevant results when pos¬ 
sible, but we make no effort on reviewing all of them 
considering the huge amount of literature. The paper 
is organized as follows: In Sec. II, we give a brief in¬ 
troduction of the EPI problem starting from the Born- 
Oppenheimer approximation. We then introduce the sys¬ 
tem setup and Hamiltonian we use in this paper. In 
Sec. Ill, we briefly summarize our use of the nonequi- 
librium Green’s function (NEGF) method to study elec¬ 
tron, phonon transport and their interaction perturba- 
tively. We consider several applications of the method. 
The first one is the effects of EPI on the thermoelectric 
transport coefficients in a single level model. The second 
one is the heat transport between electrons and phonons 
due to EPI. The use of simple models enables us to ap¬ 
proach the problems semi-analytically. The last example 
is a numerical study of the Joule heating and phonon- 
drag effect in carbon nanotubes. In Sec. IV, we consider 
the case of strong EPI using the quantum master equa¬ 
tion (QME) approach. After reviewing the earlier work, 
the same thermoelectric transport model is re-visited fo¬ 
cusing on how the strength of EPI affects the results. In 
Sec. V, we focus on the current-induced dynamics. Based 
on the Feynman-Vernon influence functional approach, 
we derive the semi-classical Langevin equation, taking 
into account the equilibrium phonon and nonequilibrium 
electron baths. The final section is our conclusion and 
remarks. 


II. BORN-OPPENHEIMER APPROXIMATION 
AND ELECTRON-PHONON INTERACTIONS 

To discuss the meaning and formulation of the elec¬ 
tron and phonon systems and their mutual interac¬ 
tion, we need to start from the Born-Oppenheimer 
approximation^’^. Consider an electron-ion system with 
a total Hamiltonian H = iJ®, where P* is the ki¬ 
netic energy operator for the ions, and U 

is electron Hamiltonian with kinetic energy of the elec¬ 
trons, P^, and potential energy U = 
which includes the Coulomb interactions among the elec¬ 
trons and ions. Since the ions are much heavier than the 
electrons, one can treat the ion kinetic energy term as a 
small perturbation with the expansion parameter^ 



where tUq is the mass of an electron and rup mass of an 
ion (assuming all have the same mass). If the ions are 
considered infinitely heavy, the ions will not move and 
the electron wavefunctions satisfy 

( 2 ) 


where x represents the set of all coordinates of the elec¬ 
trons, R the positions of all the ions, and a is the elec¬ 
tronic state quantum number. The eigen-functions and 
the eigenvalues depend on R parametrically. 

We assume an ortho normal set {(j)a} that satisfies 
Eq. (2) has been obtained. To take into account the 
effects of the ions, we consider a trial full wavefunction 
in a factored form 

T(x, R) = R)x^-cx{R) = WP), (3) 

and consider the variational solution^ of the full Hamil¬ 
tonian, min^(T|P|T), subject to the normalization 
(xlx) = 1- This variational approach is equivalent to 
omitting the off-diagonal elements (which is the Born- 
Oppenheimer approximation, see Ref. 2, App. VHI), giv¬ 
ing an equation for the ions 

+ + (</.„ |p>„) 

- — (</>«! Vfl|0„)-vdx = Px, (4) 

rup / 

where (• • •) means the x-dependence is integrated out 
but still P-dependent; Vi? is a multi-dimensional gradi¬ 
ent operator with respect to R. Since the left-hand side 
depends on the electronic quantum number a, the full 
eigen-energy E and functions also depend on a paramet¬ 
rically, e.g., we may write E/ 3 ;a- 

If we assume that the electrons are in its instantaneous 
ground state, the ions move in a potential surface gen¬ 
erated by the electrons. There are no explicit electron- 
phonon interaction (EPI) terms. To account for the EPI, 
we need to go back to the basis, Eq. (3), and consider the 
matrix elements 

(a/?|P|a'/?'). (5) 

The off-diagonal terms are interpreted as the EPI^’^, 
which are small. If the off-diagonals are omitted, the elec¬ 
trons stay in a given quantum state a. The off-diagonal 
terms describe the scattering of the electrons to different 
state a\ If ion displacements are small, the most impor¬ 
tant contribution is from the linear term in the displace¬ 
ment 

jR' , 

- {4>aX/3-,a\'^R\4>a') ' («,/3) ^ {a',/3'). 

TTLp 

(6) 

These off-diagonal matrix elements can be used, e.g., in 
a Fermi-Golden rule calculation of scattering processes. 
However, the identity (in the sense of effective Hamil¬ 
tonians) of the electrons and phonons and their mutual 
interaction are not at all clear. Although EPI plays ma¬ 
jor role in many physical processes^, such as electronic 
transport and superconductivity, its conceptual founda¬ 
tion is still not very solid. Within the Born-Oppenheimer 
scheme, it is not clear at all how to transform the origi¬ 
nal Hamiltonian H into a form of an electron system and 
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independent phonon system and their interaction unam¬ 
biguously. The problem is related to the fact that in de¬ 
riving the phonon Hamiltonian (the potential surfaces), 
the effect of electrons is already used. Thus, putting the 
electrons back amounts to double counting, see Refs. 4 
and 8 for some of the modern treatments. 

Instead of pursuing a self-consistent theory of EPI from 
the Born-Oppenheimer approximation, here in this re¬ 
view, and also in many of the practical applications^^“^^, 
we adopt a phenomenological point of view, and use the 
model Hamiltonians as given below in Eqs. (8) and (10). 
Eocusing only the term linear in the displacements away 
from the equilibrium positions of the ions, we can think 
of the single electron Hamiltonian below having a R- 
dependence. Taylor expanding it, R = R? u /, we 
obtain 





' ' dRk 


( 7 ) 


where \j) is the single particle state when ionic sys¬ 
tem is in equilibrium position R^. The extra factor of 
square root of ion mass is because of our conven¬ 
tion of displacement variable u. This form of interaction 
is intuitively understandable and originally proposed by 
Bloch^^. In Chap. 4 of Ref. 6, a derivative from Eq. (6) 
to (7) is given, but the reasoning does not seem to be 
rigorous. 

Thus, our starting point of a derivation is a tight- 
binding Hamiltonian for the electrons, harmonic cou¬ 
plings for the phonons and a standard EPI term. They 
are taken as given and exact. The charge redistributions 
and self-consistency for the electrons are not part of the 
discussion. Symbolically, the total many-body Hamilto¬ 
nian is given as 


( 8 ) 


where the electron part is = c^Hc^ the phonon 
part i7p = + u^Ku) + Vn{u^)‘ The variable 

u is mass normalized, Uj = yRn^{Rj — R^)- Because 
of this, the conjugate momentum is p = ii. is 

the nonlinear force contribution, c is a column vec¬ 
tor of the electron annihilation operators, which we can 
separate into three regions, the left, center, and right, 
c = (c^,d, c^)^, T stands for matrix transpose. Sim¬ 
ilarly u = {u^, . Accordingly, the matrices H 

and K are partitioned into nine regions (submatrices), 
e.g., 


/ 0 \ 

H = , (9) 

\ 0 ) 


such that ^ -h Re, with K = + 

Vj" = + H.C.. Note that we assume no 

interaction between the left and right leads (See Ref. 67 
for transport when there is a lead-lead coupling). We do 


a similar partitioning for K using the notation of Ref. 68. 
The EPI takes the form 

ifep(rf, U^) = = E (10) 

ijk k 

We assume that the EPI appears only in the central re¬ 
gion. A schematic representation of the system setup is 
shown in Eig. 1. 

The separation of the electron and phonon leads makes 
the theoretical development easier. In reality, they could 
either be physically separated, or built into one. Eor ex¬ 
ample, one electrode could serve both as an electron and 
a phonon lead, but we assume that we have independent 
control over their temperatures and T^, a = L^R. 



FIG. 1. Model system considered in this review. The center 
device, including both electrons, phonons, and their interac¬ 
tions, is coupled with two electron and two phonon leads. 
Each electron lead is characterized by its chemical potential 
fia and temperature ^ and each phonon lead by tempera¬ 
ture Tp^. 


III. WEAK EPI REGIME: NONEQUILIBRIUM 
GREEN’S FUNCTION METHOD 

A. Theory 

We first consider the case where EPI is weak, so that 
we can perform a perturbation expansion over the inter¬ 
action matrix M. In order to do so, we use the NEGE 
method. Detailed introduction is given in our previous 
work^^’^^“^^. This section can be considered as an appli¬ 
cation of the general approach developed in Refs. 68-69 
to the EPI problem. We use similar notations therein, 
and only give a brief outline of the approach here. 

We denote the electron device Green’s function with¬ 
out and with EPI by Go and G, the corresponding 
phonon Green’s functions by Dq and D, and the lead 
Green’s functions without coupling to the center as pa 
and dc, respectively. The couplings of the device with 
the leads and that between the electrons and phonons 
are described by self-energies, with E and H representing 
that of electron and phonon, respectively. Eor example. 
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we define the time-ordered electron Green’s function in¬ 
cluding EPI on the Keldysh contour [Fig. 14 (b)] 

G'ii(r,T') = -^(rcCi(T)ct(r')). (11) 

Here, r/r' is time on the contour, and i/j is index of the 
electronic states. The contour time order operator Tc 
puts the operators later in the contour to the left. The 
average (•) is with respect to the density matrix of the 
full Hamiltonian. 

The contour ordered Green’s function can be divided 
into different groups according to the spatial position of 
i/j, similar to the Hamiltonian. The most interesting 
one is G^, where i and j are both at the center device 
region. At the same time, it can be written as a 2 x 2 
matrix in time space 


the interaction Hamiltonian i^ep, using Feynman diagra- 
matics. The interacting Green’s functions are expressed 
using similar Dyson equations as Eqs. (13-14), 

C?(l,2) =Go(l,2) + Go(l,3)I]ep(3,4)G(4,2), (15) 
D(l, 2) = Z)o(l, 2) + Do(l, 3)nep(3,4)D(4,2). (16) 

Here, Fep and Hep are electron and phonon self-energies 
due to EPI. Using Eq. (12), at steady state, we can get 
the following useful relations in energy/frequency domain 

G^(£) = [(£ + i0+)/ I]Lt(e)] , (17) 

D'-iuj) = [(u; + iQ+fl -K^- n[,t(a;)] , (18) 

S[,t(£) = S[(£)+I];p(£), (19) 

nL,(a;)=nUw)+n;p(a;). (20) 


G{Ti,Tj) 




( 12 ) 


with G^, G^, G^, G^ the time-ordered, anti-time- 
ordered, greater and lesser Green’s functions. The re¬ 
tarded and advanced Green’s functions are obtained from 
them, i.e., G^ = G^ - G<, and G^ = G< - GK For the 
definition and relations among these Green’s functions, 
we refer to the book by Hang and Jauho^^, and our pre¬ 
vious publications^^’^^’^^. 

To calculate the Green’s functions, we use a process of 
two-step adiabatic switch on. We start from the decou¬ 
pled system and leads. Each of the electron and phonon 
leads is at its own equilibrium state, characterized by 
the temperature T" and/or chemical potential fia- The 
corresponding equilibrium Green’s functions can thus be 
defined according to the equilibrium canonical distribu¬ 
tion. The initial state of the system is arbitrary and not 
important in most cases (e.g., for steady state). 

At the first step, we switch on the interaction of the 
center Hamiltonian with the electron and phonon leads. 
We wait until the electron and phonon subsystem reaches 
their own nonequilibrium steady state, since the temper¬ 
ature and/or chemical potential of each lead can be dif¬ 
ferent. The two subsystems are quadratic and exactly 
solvable, and we get the non-interacting center Green’s 
functions Gq and Dq from the Dyson equation (we omit 
the superscript G) 

Go(l, 2) = gcih 2) + gcih 3)Sb(3,4)Go(4,2), (13) 
Do(l, 2) = dc(l, 2) + dcil, 3)nb(3,4)Do(4,2). (14) 

Here, we have used a single number to represent the ma¬ 
trix indices and contour time arguments, i.e., Go(l,2) = 
Summation or integration over repeated 
indices is assumed. gc {dc) is the center electron 
(phonon) Green’s function without coupling to the L and 
R leads. The self-energy Eb = Ei, -1- E^^ includes contri¬ 
butions from L and R, with Eq,( 1, 2) = 
similarly for Hb- 

At the second step, we adiabatically switch on the EPI 
in the center. We perform a perturbation expansion over 


We use 5 for the energy of electron and uj for the angular 
frequency of phonon, respectively, and I is the identity 
matrix. To get an expression for the current, we also need 
the greater and lesser version of the Green’s functions^^ 

G>’<(£) = G'-(£)S>’<(£)G“(£), (21) 

D>’<{co) = D^{co)Il>^<ico)D‘^ico). (22) 


The electrical current (R) is expressed as the change rate 
of the electron number in one of the leads times the 

charge of electron (—e). For example. 




dt 

= — — ImTr 
a 


{c^\t)d{t)) 


= 2eReTr = 0)] . 


(23) 


It can be expressed by the Green’s function of the center 
region and the lead self-energies^^“^^. 



(24) 


Similarly for the heat current carried by electrons (Jh) 
and phonons (Ip) 

h = \ liL)Tr [G>S< - G<S>] , (25) 

/ + 00 7 

^ [D>n< - Z)<n>] . (26) 


We have defined the positive current direction as elec¬ 
trons going from the lead to the center. We dropped 
the argument of the Green’s functions for simplicity. We 
ignore the spin degrees of freedom, since it is not rele¬ 
vant here. Gurrents out of the right lead are obtained by 
replacing index L by R. One can symmetrize the expres¬ 
sions based on energy and charge conservation. 

The set of coupled equations Eqs. (17-22) is difficult 
to solve, due to the many-body EPI. Since the EPI is 
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weak, we consider only the lowest order Feynman dia¬ 
grams shown in Fig. 2. The expressions for the self¬ 
energies are as follows. The electron Fock self-energies 
from phonons are 

(27) 

+ G>,/£_)DoXa;))My^. (28) 

The Hartree self-energy does not depend on energy 

SLn = = 0) I(29) 

The phonon self-energies from electrons are 

n<’>(a;) = [M'"G<’>(e)M"G>’<(£_)] , 

(30) 

nrn(^) = -*/ [M™Gy(e)M”Go<(£-) 

+ M'"G<(e)M"G“’''(e_)] , (31) 

with 5- = e — hjj. Summation over repeated indices is as¬ 
sumed here. Different charge and energy conserving ap¬ 
proximations have been developed in the literature. We 
will use two of them. In Subsec. IIIB, we perform an 
expansion of the current up to the second order in M, 
following the idea of Ref. 29. In the numerical model cal¬ 
culation in Subsec. IIIC, we use the self-consistent Born 
approximation (SCBA), which means we replace Go, Dq 
by G, D respectively in the above equations. 



FIG. 2. Feynman diagrams due to electron-phonon interac¬ 
tion. The first two are Hartree and Fock diagram for electrons, 
and the last one is the polarization bubble for phonons. The 
expressions of these diagrams can be found in Eqs. (27-31). 


B. Thermoelectric transport through a single 
electronic level 

We consider a single electronic level = Sod^d^ cou¬ 
pled to the left and right electrodes, characterized by the 


constant level-width broadening Fq, with energy cutoff 
Ed (see Eq. (40) for the general definition). It interacts 
with an isolated phonon mode with frequency cjo, and 
^ep = mod^du. In the linear regime, we introduce an 
infinitesimal change of the chemical potential or temper¬ 
ature at lead L, e.g., = /i + = T + with 

(7 = e or p, /i and T are the corresponding equilibrium 
values. We look at the response of the charge and heat 
current due to this small perturbation. The result, up to 
the 2nd order in M, is summarized as follows. 


4e 


4 


£o 

Cl C2 



(32) 


The linear conductance and the Seebeck coefficient are 



Ge — 

(33) 


(5R Cl 

(34) 




5T~ eCoT' 

The coefficients Cn are 



3 

i=l 

(35) 

with 



/:(i) - 


(36) 

ri^) - 


(37) 

£(3) _ 

J 1 J(e - m)” (A7i;(r„ + 2GSImE^pGgr„ 

)/'• 



(38) 

We have defined 



II 

1 

1^ 

(39) 


r„ = z(i];-s“), 

(40) 


A /^r-pe /^a 

— 'Argi 

(41) 


= GSReE^pAa + 7laReE“pGg, 

(42) 


AA" = zGSImE^pAg + zAgIniE“pGg, 

(43) 


and a means the lead different from a. Cn^ is the single 

('2') 

electron Landauer result. Cn ^ is the quasi-elastic term. 
£(3) 

is the inelastic term. / is the Fermi-Dirac distribu¬ 
tion function 


/a(e) 


exp 


/ g - /ia \ 

V ) 



(44) 


Since we are looking at the linear response regime, /l = 
/i^, we dropped the subscript in Eq. (39). We will also 
use the Bose-Einstein distribution later 


«b(w,T) 



(45) 


When there is no ambiguity, we will also drop the argu¬ 
ment T. 
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In the following, we set the position of the electronic 
level to 5o = 0, and look at the dependence of the con¬ 
ductance on the chemical potential ja. We firstly write 
down the expressions for the self-energies, and make some 
observations based on their functional forms. 

The Hartree self-energy is real and does not depend on 
energy 


^ 


27rujn 


L 


0 J —Ed 


faje) 

2 + rV4 


de. 


(46) 


At T = 0, we get 


^ 


7TUJ. 


0^ 


_.2e 
tan — 


(47) 


with r = Tl -\-rji. For large enough sd^ the lower limit 
term turns to —mo/(2co’o), which is the polaron energy 
shift. Note here that the 1/2 is due to the fact we use 
in our definition of (Eq. 10). This is different 
from the common definition that uses the creation and 
annihilation operatore a’*' + a (Eq. 74). We have sub¬ 
tracted the polaron shift term in the following calcula¬ 
tion, since it is a constant. After this subtraction, 
is odd in fi with a negative slope near /i = 0. We fo¬ 
cus on the /i > 0 regime. It saturates to —mp/{2ujp) for 
large e.g., At non-zero temperature, the slope 

and the saturation value change, but the shape of the 
curve is similar to the T = 0 case. This means that, the 
Hartree term shifts the electronic level, and reduces the 
conductance. On the other hand, when /i ^ T, the con¬ 
ductance tends to zero, whether we include the EPI or 
not. Thus, the correction to the conductance due to the 
Hartree term AGeH <0- It starts from zero at /i = 0, 
goes back to zero at /i ^ T, and it reaches a maximum 
magnitude at some point in the middle. The described 
behaviour is schematically shown in Eig. 3 (a). This effec¬ 
tively reduces the broadening of the single level spectral 
function, also the conductance peak. Since the Seebeck 
coefficient is related to the logarithmic derivative of the 
conductance, we expect it to increase the magnitude of 
the Seebeck coefficient near resonance, and to reduce it 
off resonance [Eig. 4]. 

The imaginary part of the retarded Eock self-energy 
is32 

2 

Iml]^(e) =W sAa{e - sfkoo) 

x[l+nBisu)o) - fais - shvo)]. (48) 

It is negative and even in 5. Its role on the differential 
conductance at the phonon threshold {eV = hujo) has 
been discussed extensively^^“^^. The main conclusions 
are: it reduces the differential conductance at eV = hwo 
for resonant case (/i ^ 0), where the bare transmission 
without EPI (To ^1), while it does the opposite for far 
off resonance case (/i ^ F), where Tq ^ 0. The transition 
point between the two opposite behaviors is Tq = 1/2 if 
the electronic density of state (DOS) is fiat. But, in 


general, it depends on the system parameters. At non¬ 
zero temperature, the sharp threshold broadens out, and 
the linear conductance is affected: its correction to the 
conductance AGeFi is negative for /i ^ 0, and positive 
for /r > F [Fig. 3 (c)]. 

Physically, ImE^ gives rise to phonon scattering pro¬ 
cesses. Its effect can be understood as follows: At T = 0, 
for small bias {eV < huJo), phonon emission is not possi¬ 
ble due to Pauli blocking, while phonon adsorption is not 
possible due to zero phonon population. So, ImE^ does 
not affect the linear conductance. At high enough tem¬ 
perature, both phonon emission and adsorption are possi¬ 
ble even at small bias, due to the broadening of the Fermi 
distribution, and finite population of phonon modes. The 
phonon scattering process decreases the conductance on 
resonance, but increases it far off resonance. As a result, 
the Seebeck coefficient becomes smaller. 

The real part ReE^(5) is obtained by the Hilbert trans¬ 
form of the imaginary part. At zero temperature, it di¬ 
verges logarithmically at 5 — /i = . Its effects on 

the conductance and Seebeck coefficient are difficult to 
analyze. We rely on the numerical result [Fig. 3 (b)]. 

Figure 3 shows the correction to the linear conductance 
of different self-energy terms as a function of ja. These 
numerical results confirm our qualitative analysis. By 
comparing the total conductance at low [Fig. 3 (e)] and 
high temperature [Fig. 3 (f)], we see that, (1) the Hartree 
term dominates at low temperature, and the Ge-/^ peak 
becomes narrower. (2) the Fock term becomes important 
at high temperature, and the Ge-/i peak broadens out. 
Their effects on the Seebeck coefficient (S) are shown in 
Fig. 4. At low temperature, when the EPI is included, 
the magnitude of S gets larger for /i ^ 0, and smaller 
for \fi\ ^ F. At high temperature, the effect of the Fock 
term results in drop of S. In any case, the correction to 
S is small for weak EPI. But for the case of strong EPI, 
the correction could be large (see Subsec. IV C). 


C. Heat transport between electrons and phonons 

Let us look back at the setup in Eig. 1. We want to 
study the heat transport between electrons and phonons 
at finite temperature bias, but zero voltage bias. The 
simplest setup is that the system couples to one electron 
and one phonon lead, each at its own temperature, see 
Eig. 5. The expression for the energy current from elec¬ 
trons to phonons can be obtained from Eq. (26) and the 
expressions of the self-energies Eqs. (27-31)^^ 

^ = [G>„(£)M^,Z)<(a;)G<.(£_)Mj„] , 

(49) 

where, again, summation over repeated indices is as¬ 
sumed. Eor the ease of analysis, we perform an expansion 
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FIG. 3. Chemical potential dependence of the conductance 
correction due to the ReE5^(5), ImE5n(c), and the sum 

of them (d). The blue, red, purple, and black lines correspond 
to ksT — 0.05,0.15,0.25,0.35da;o, respectively, (e)-(f) The 
electrical conductance as a function of chemical potential (/x) 
at ksT — 0.05dct;o (e) and ksT — 0.35da;o (f). Solid lines 
include EPI, while the dashed lines do not. Other parameters 
used: Tl = Ti^ = — 1, mo = lda;o/(AyT[). 
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FIG. 4. Change of the Seebeck coefficient due to EPI at dif¬ 
ferent temperatures. The parameters and meaning of colors 
are the same as Fig. 3. 


of the above expression to 2 nd order in M, and it becomes 



FIG. 5. The model system we consider to study the energy 
transport between electrons and phonons. 


where 

A(w,Te) = j ^Tx[MA{e)MA{e-)\ 

X [/(£,Te) -/{£_, Te)], (51) 

and 

A{u:)=i[Dl{e)-D^,{e)\. (52) 

Now the question we ask is whether there is a diode 
behaviour for the heat transport between electrons and 
phonons^^, e.g., Q{AT) ^ Q{—AT)^ with AT = — Tp. 

This is relevant because in some special situation, e.g., 
at metal-insulator interface, or insulating molecular junc¬ 
tions, EPI becomes the bottleneck of heat transfer. 
We can define the rectification ratio as 

_ Q{AT) + Q{-AT) 

Q{AT) - Q(-AT) • ^ ’ 

If we assume a constant electron DOS, A(ci;) does not 
depend on T, we get Q{—AT) = —Q{AT), and R = 0. 
The physical reason is that in this case, it is possible 
to map the electron-hole pair excitation into harmonic 
oscillators^^’^^’^"^. Then, it is equivalent to heat trans¬ 
port within a two-terminal harmonic system. We do not 
expect any rectification effect. To make i? 7 ^ 0, the elec¬ 
tronic DOS has to be energy-dependent within the broad¬ 
ening of the Fermi-Dirac distribution given by ksAT. 
This effectively introduces anharmonicity into the sys¬ 
tem, consistent with previous studies^^“^^. 

We can go one step further, by making a Taylor ex¬ 
pansion of the spectral function A{s) about the Fermi 
energy, we find that the sign of R is determined by the 
sign of (The 1st order term is zero). 

To check this argument, we calculated the heat cur¬ 
rent across one-dimensional (ID) metal-insulator junc¬ 
tion. The metal side is represented by a ID tight-binding 
chain, with hopping element t = —0.1 eV, and onsite 
energy £0 = 0- The insulator side is represented by a 
ID harmonic chain with the spring constant K = 0.1 
eV/(A^u). The insulator and metal couple through their 
last two degrees of freedom. Their interaction matrices 
are 


r-\-oo j 

^ [A(w, Te).4(w)] 

Jo 

X [nB(w,re) - nB(w,rp)], 


= (-l)'=mo 5 J ) , fc = l,2. (54) 

Here, mo = 0.05 eV/(AA/u), k represents the phononic 
(50) degrees of freedom. This means that the system couples 
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FIG. 6. (a) The energy current with /i = 0 (red, solid) and 
/i = 0.1 eV (blue, dotted). We have set — T -\- /ST 12^ Tp = 
T — AT/2, with T = 300 K. (b) Fermi level (/x) dependence 
of the energy current at fixed temperature difference AT = 
300 K. (c) Rectihcation ratio T as a function of /x. (d) The 
electronic spectral function A{e) as a function of energy e. 


to only one electron and one phonon lead [Fig. 5]. Fig¬ 
ure 6 summarizes our result. In Fig. 6 (a), we show the 
energy current as a function of temperature difference be¬ 
tween the metal and insulator (AT), at two Fermi levels 
(/X = 0,0.1 eV). The energy current is asymmetric with 
respect to the sign change of AT. But the rectification 
ratio R has opposite sign. This is further highlighted in 
the plot of the energy current and the rectification ratio R 
as a function of /x for fixed temperature bias AT = ±300 
K [Fig. 6 (b) and (c), respectively]. The sign of R is well 
correlated with the sign of [Fig. 6 (d)]. 


D. Effect of EPI on thermal transport in 
single-walled carbon nanotubes 

Phonon modes can be excited by the mobile electrons 
due to the EPI effect; i.e., a high bias over the system 
leads to self-heating (Joule heating). In nanoscale electric 
devices, the electric current density can be much larger 
than that in the macroscopic system. The high current 
density will generate strong Joule heating, which may 
eventually break the device. In this sense. Joule heating 
becomes a bottleneck for further increase of the electric 
current density. Hence, lots of theoretical and experimen¬ 
tal efforts have been devoted to understanding the Joule 
heating phenomenon in the nanoscale electric devices. In 
experiment. Joule heat can be measured via the thermal- 


mechanical expansion technique, which records the Joule 
heat induced temperature rise^^. The Joule heat can in¬ 
crease the temperature of the molecular junction from 
room temperature to 463 K, which has been examined 
through the inelastic electron tunneling spectroscopy^^. 
Grosse et al. investigated the nanoscale Joule heating 
in phase change memory devices^^. The Joule heating 
leads to the temperature rise in the phase change mem¬ 
ory device, which results in an obvious volume expan¬ 
sion. In another experiment. Joule heating is found to 
be responsible for the correlated breakdown of nanotube 
forests^^’^^. 

For a system without localized phonon modes, all 
phonon modes have important contribution to the Joule 
heat. The Joule heat contributed by these propagating 
phonon modes have important effects on the electric de¬ 
vices. For example, in graphene transistors, the output 
electric current will saturate with increasing source-drain 
voltage^^’^^. This saturated current density can be re¬ 
duced by 16.5% due to the Joule heating^^. 

The localized phonon modes exist around some defects 
or nonuniform configurations, such as the free edge, the 
isotropic doping, interface, etc. This particular type of 
phonon modes has no direct contribution to the thermal 
conduction, but localized phonon modes play a particu¬ 
larly important role in the Joule heat phenomenon. They 
are characteristic for their exponentially decaying vibra¬ 
tion displacement; i.e., only a small portion of atoms 
are involved in the localized vibration. For instance, 
there are some localized edge phonon modes at graphene 
nanoribbon’s free edge. In these modes, edge atoms vi¬ 
brate with large amplitude, but the vibrational displace¬ 
ment decays exponentially from the edge into the interior 
region. 

The localized-phonon-mode-induced Joule heat was 
observed in graphene nanoribbons in experiment, and 
explained theoretically. Jia, et al. utilized Joule heat¬ 
ing to trigger the edge reconstruction at the free edges in 
the graphene nanoribbons^^. Engelund, et al. attributed 
this phenomenon to the Joule heating of the edge phonon 
modes^^. There are two conditions for the important 
Joule heating of the edge phonon modes. Eirst, these lo¬ 
calized edge phonon modes can spatially confine the en¬ 
ergy at the edges. Second, the electrons interact strongly 
with the localized edge phonon modes. The mean steady- 
state occupation of the edge phonon mode can be calcu¬ 
lated from the ratio of the current-induced phonon emis¬ 
sion rate and damping rate. The effective temperature 
for the free edge can be extracted by assuming this occu¬ 
pation to be Bose distributed. The effective temperature 
was found to be as high as 2500 K for bias around 0.55 V. 
This high effective temperature was proposed to be the 
origin for the edge reconstruction. 

Although Joule heating might be used for selectively 
bond-breaking^^its most common outcome is a disas¬ 
ter of device breakdown. The effective temperature is a 
suitable quantity to describe the Joule heating. In 1998, 
Todorov studied Joule heating problem in a molecular 
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FIG. 7. (Color online) The phonon thermal conductance ver¬ 
sus chemical potential for metallic SWCNT (10, 10) at (a) 
150 K and (b) 300 K. Solid line is for ballistic phonon ther¬ 
mal conductance without EPI effect. Reprinted from J. Appl. 
Phys., 110, 124319 (2011). 


junction^^. In his work, the Einstein model is applied 
to represent the phonon modes in the system, and the 
electron-electron interaction is ignored as the system size 
is much smaller than the electron mean free path. Part 
of the EPI-induced Joule heat will be delivered out of 
the system by the phonon heat conduction, while the re¬ 
maining Joule heat gives a high effective temperature. 
For low ambient temperature, the effective temperature 
scales with voltage V as ^ with 7 as an EPI- 

dependent constant. It was shown that the effective tem¬ 
perature can be above 200 K for a very low ambient tem¬ 
perature around 4 But at very high bias, the scaling 
law could differ from this^^. 

There are several experimental approaches to inves¬ 
tigate the effective temperature of the electric device 
induced by Joule heating. The effective temperature 
can be extracted by measuring some quantities that are 
temperature-dependent. For example, the breaking force 
of the single molecular junction is related to the tempera¬ 
ture. This force-temperature relationship can be used to 
estimate the effective temperature^^. The Raman spec¬ 
troscopy also depends on the temperature. Hence, it can 
be used to deduce the effective temperature of Raman- 
active phonon modes^^’^^’^^. 

It has also been shown that EPI has an important 



p(eV) 


EIG. 8. (Golor online) The phonon thermal conductance ver¬ 
sus chemical potential for semiconductor SWGNT (10, 0) at 
(a) 150 K, (b) 300 K, and (c) 1000 K. Solid line is for ballistic 
phonon thermal conductance without EPI effect. Reprinted 
from J. Appl. Phys., 110, 124319 (2011). 


effect on the thermal conductance in single-walled car¬ 
bon nanotubes (SWCNTs)^^. For them, we apply the 
Born approximation to consider the EPI effect using the 
NEGF approach, as the SCBA is computationally more 
expensive. The phonon thermal current can be calculated 
by considering the three EPI contributions shown in the 
Feynman diagrams in Fig. 2. The phonon thermal cur¬ 
rent flowing from the left lead into the center is given by 
Eq. (26). The expression for the right lead is analogous. 
The Joule heat is generated in the system and flows into 
the leads, so the total Joule heat is the sum of heat cur¬ 
rents into both leads, Q = — (Jp +^^)- The thermal cur¬ 
rent from Eq. (26) also includes that induced by the tem¬ 
perature gradient, which satisfies Hence, Q 

gives solely the Joule heat.For metallic SWCNT (10, 10), 
both electrons and phonons are important heat carriers. 
The EPI only slightly reduces the electron thermal con¬ 
ductance, but it has a strong effect on the phonon ther¬ 
mal conductance. More specifically. Fig. 7 (a) shows an 
‘electron-drag’ effect on the phonon thermal conductance 
at 150 K. The phonon thermal conductance becomes neg¬ 
ative for high chemical potential value /i > 2.0 eV, which 
indicates that electrons can help to drag phonons from 
cold temperature region to the hot temperature region. 
The ‘electron-drag’ phenomenon happens at low temper¬ 
ature and high chemical potential, and it does not hap¬ 
pen at a higher temperature 300 K as shown in Fig. 7 (b). 
For semiconductor SWCNT (10, 0), the electronic ther¬ 
mal conductance contributes less than 10% of the to¬ 
tal thermal conductance at low bias (e.g. /i = 0.3 eV), 
while phonons make most significant contribution to the 
total thermal conductance. Similar ‘electron-drag’ phe- 
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nomenon also exists in the semiconductor SWCNT (10, 
0) at low temperature as shown in Fig. 8 (a). 


IV. STRONG EPI REGIME: QUANTUM 
MASTER EQUATION APPROACH 


A. Quantum Master equation formulism 

In this section, we introduce the QME approach to 
consider the case of strong FPL Before doing that, we 
should mention that the NEGF method has also been 
used to treat the strong Since the idea behind 

it is very similar to that of the master equation approach, 
we choose not to introduce it here. 

To simplify the formula, we ignore the coupling of 
molecular phonon modes to the phonon leads. The model 
Hamiltonian simplifies to 

(55) 

where Hs = denotes the system Hamil¬ 
tonian, and Ve = is the system-lead coupling. 

In the QME formalism we assume the system-lead cou¬ 
pling Ve is weak so we can do perturbation on it. We 
work in the interaction picture with Hq = H — Ve diS 
non-interacting part and Ve as the interaction. Eor sim¬ 
plicity, in this section we use V to represent Ve since we 
don’t have Vp. The equation of motion for the full density 
matrix follows the von Neumann equation 


dpijt) 

dt 




(59) 

where C^^{t — t') = Tr[pg (g) is the cor¬ 

relation function of the leads. Here we have used the 
condition that the expectation value of a single lead op¬ 
erator B^ is zero. We can now transform back to the 
Schrodinger picture and extend the initial time to to — oc 
to get the QME of Redfield type^^^“^^^ 


a,/3 

X f dt'[[S^, {t' — t)p]C"^(t — t') H- H.c.}. 

J —oo 


Here we have replaced p(to) by p, which is essential 
and correct only when one intends to get the 0th- 
order reduced density matrix by solving the above 
In the application to the EPI problem, 
by exact diagonalizing the system Hamiltonian, this Red- 
field QME can take into account the coherence between 
electrons and phonons, in contrast to the usual rate equa¬ 
tion approach^^^’^^^. 

We write the above equation in the eigenbasis of the 
system Hamiltonian Hs lo obtain^^^’^^^ 

= + (61) 


ih—^ = [Vi{t),pi{t)]. (56) 

Here, the subscript / denotes operator in the interaction 
picture. The time argument in the parentheses means 
non-interacting evolution 0{t) = Xhe 

above equation can be written in an integral form as 

Pi{t) = dt'[Vi{t'), Pi{t')] + Piito). (57) 

One can recursively apply the above equation to get a se¬ 
ries expansion of the full density matrix in power of Vj. 
We truncate the series to the second order and differenti¬ 
ate it with respect to time t at both sides of the equation 
to get the following integro-differential equation 

^^y[n(i),P/(io)]-^^‘rfi'[n(i),[n(iO,/^/(io)]]. 

(58) 

We prepare the initial state as a product state of the 
system and each lead, p{to) = p{to) ^ Pe ^ P^‘ For the 
system-lead coupling V we assume it can be written as 
a product of system operator S and lead operator B as 
V = S'^ (g) B ^. In such cases we can trace over the 

lead degrees of freedom to get 


where the relaxation tensor reads^^^ 
mj _ Jl qd 

a,(3 

(62) 

I 

The transition coefficients are given by 

W^f = f (63) 

J —OO 

where Akj = — Ej are the energy spacings of the 

system Hamiltonian. 

Since we are only interested in the steady state, we 
impose the condition dp/dt = 0 at t = 0 and solve the 
above equation order by order with respect to V. One 
can find that all the off-diagonal elements of the 0th- 
order reduced density matrix vanish in steady state and 
the diagonal elements can be evaluated via the matrix 
equation^^^ together with the constraint 

of Tr[p^o)] = 1. 

Eor the calculation of currents, we go through a sim¬ 
ilar derivation as the QME. The electronic current op¬ 
erator Je and heat current operator can be writ¬ 
ten in the form J'e(h) = ^ '®e(h)' "Fhe expec¬ 

tation value of currents can be calculated according to 
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^e(h) = Tr[p/(t) Since we are interested in the low¬ 

est order of current, we truncate Eq. (57) to the lowest 
order and plug in to get 

4(h) = Y / dt'TT{[Vi{t'),piito)]Jiit)}. (64) 


By taking the trace over the leads and transforming back 
to the Schrodinger picture one obtains the current at t = 
0 as 


4(h) = ^ 


a,P d-oo 


) + H.c., 


(65) 


where = (-B“(i)'Bg^h)(*^)) correlation func- 

tion between the lead operators occurring in the system- 
lead coupling Hamiltonian and the current operator def¬ 
inition. For the same reason as the derivation of master 
equation, here p{to) needs to be replaced by p to get cor¬ 
rect steady state results. Since we are calculating lowest 
order of current, we can use the 0th order reduced density 
matrix p ^^^. Written in the eigenbasis of system Hamilto¬ 


nian, the above equation becomes /e(h) = Tr 
with the reduced current operator defined as 

(4V))^i = p E [5a4>Ve“(h)(Afe4+c.c. 




h? ^ - 

a,l3,k 


'e(h) 


( 66 ) 


where the transition coefficients are 


KU^kj) = f (67) 

J —OO 


Up to now the QME formalism is general and not 
restricted to any specific form of system or leads 
Hamiltonians. For the application to transport prob¬ 
lems with EPI as concerned in this review, the 
system-coupling Hamiltonian is considered as a tun¬ 
neling Hamiltonian^^. In such case the system op¬ 
erator S', leads operator B and B can be specified 

as S = {d, d”*"}, B = |Z]/cGL,i? Z]/ceL,i? 

= (eEfegLUct-eEfeeL^feCfc I and Bh = 

{ - T,keLi^k - fJ-LWkCk |. The in- 

finite nature of the leads can be specified by defining a 
continuous spectra function for the leads 

r„(£) = -2ImE;(£) 

= 2wJ2\Vk\'^H£-£k), a = L,R. (68) 

k£a 

Throughout this section we use a wide band spectra func¬ 
tion for the electronic leads with Lorentzian cut-off as 


ra(e) 


Va 

l + (e/£D)2’ 


o, — T, R. 


(69) 


The non-vanishing correlation functions can be evaluated 
via 

/ OO 7 

^^^„(£)^(£)e-*/^ (70) 

/ OO 7 

^^r„(£)(l-/„(£))e--‘/'( (71) 

/ OO J 

— (e-ML)rL(e)/L(e)e‘®‘/'^, (72) 

/ OO 7 

^^4-ML)rL(e)(l-/L(£))e--*^ 

(73) 

and C^Ht) = ClHt) = -eCl^t), 

= eCj}{t). We note that here the upper index 
1 or 2 refers to the two components of B and Se/h given 
above. For the system Hamiltonian, we focus only on the 
single electronic level coupled to a single phonon mode 
representing the center of mass of the molecule. In such 
case, the system Hamiltonian will reduce to 

Hs = Sod^d + hujQa^a + Ad^d(a^ + a) (74) 

with ujQ the angular frequency of the phonon mode and 
A denotes the EPI strength. This type of Hamilto¬ 
nian has been well-studied in the context of molecular 

The QME formalism treats the nonlinearity of EPI ex¬ 
actly. Therefore in this section we will focus on strong 
EPI regime with emphasis on the EPI strength depen¬ 
dence of the electron/heat current. In the following we 
will discuss the effect of EPI on the electronic transport 
properties, including the phonon sidebands, negative dif¬ 
ferential resistance, thermoelectric properties and local 
heating effects. 


B. Phonon sidebands and negative differential 
resistance 

One of the earliest findings of the vibrational effects 
on the electronic transport through a molecular quan¬ 
tum dot is the appearance of the phonon sidebands in 
the I — V characteristics. When electrons transport 
through a single electronic level, the differential conduc¬ 
tance (dl/dV) will manifest a peak at resonant level 
when plotted against the voltage bias (U). However, 
when the electronic level is interacting with a vibrational 
mode, replica side peaks will appear at the side of the 
resonant peaks. These side peaks are called phonon side¬ 
bands. A simple reason for the appearance of phonon 
sidebands is due to the fact that the electrons can emit 
or absorb phonons when they pass through the molecule. 
Therefore, the distance between each adjacent peaks is 
always equal to a single phonon energy. The phonon 
sidebands attract wide attentions in molecular junction 
systems. Experimentally, phonon sidebands were found 
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in 1980s^^^ and then were utilized to identify vibra¬ 
tional modes in molecular junctions^^^”^^^ and quantum 
wires^^’^^“^^’^^^. Theoretically, at the beginning the side¬ 
bands were investigated by using scattering theory, which 
gives the transmission probability T(e^e') for an electron 
to passing through an EPI system. The electron is com¬ 
ing from vacuum at energy 6 and leaving at energy 
The scattering theory predicts side peaks in the transmis¬ 
sion probability, which qualitatively justifies the phonon 
sidebands in molecular junction systems. However, pre¬ 
diction of phonon sidebands in the lead-molecule-lead 
junctions is a much more difficult task. A simple gen¬ 
eralization to take the Fermi-Dirac statistics nature of 
the electron leads into account is to weight the exact 
transmission probability with the Fermi-Dirac distribu¬ 
tion of each lead, i.e., by multiplying the transmission 
probability T{e,e') by a factor fL{e)[l — fnie')] as a 
new transmission probability for electron going from the 
left lead to the right lead through the nano-conductor. 
This approximation is called single particle approxima¬ 
tion (SPA). Plenty of earlier work is in this framework 
and it predicts Lorentzian type of phonon sidebands with 
the same width. But obviously this method assumes each 
electron transports independently through the junction, 
where the many-body effects are ignored. As a result, it 
overestimates the currents and it is not able to predict 
the quantized conductance jh either. 

Based on the NEGF technique, more rigorous meth¬ 
ods merged in dealing with the nonlinearity in EPI, 
such as the Green’s function equation of motion method 
(EOM)^^, and nearest neighbor crossing 

approximation (NNCA)^^^. All of the above are Green’s 
function based formalisms with different kinds of approx¬ 
imations. In general these approaches predict that the 
phonon side peaks are much sharper than the SPA ap¬ 
proach. This sharpness is closely related to the Pauli 
exclusion, which the SPA approach failed to take into 
account. Other than Green’s function based methods, 
another approach is to use rate equation of electron oc¬ 
cupation probability in the molecule, via calculating the 
transition coefficients of the electron to tunnel from the 
molecule to each lead and vice versa. This method as¬ 
sumes the transport is an electron tunneling process and 
the electron will lose its phase information when it re¬ 
sides in the molecule. Therefore it will be valid when 
the molecule-lead coupling is weak and the coherence 
of the electron and phonon in the molecule can be ne¬ 
glected. For all these formalisms, we would like to point 
out that one should take care of the phonon distribu¬ 
tion. Treating the phonon at equilibrium distribution at 
a fixed temperature could be valid when the EPI strength 
is much weaker than the coupling strength between the 
phonon and its environment. However, when the environ¬ 
mental influence is weak, one should consider phonons in 
nonequilibrium states. This nonequilibrium treatment of 
phonon distribution can have pronounced effect on I — V 
characteristics due to the fact that the current induced 
vibrational excitation can be significant 
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FIG. 9. The phonon sidebands with different EPI strength 
with (a) low temperature (T = 0.02/ia;o/^n) and symmetri¬ 
cally biased voltage (Vl/r = zLVbias), (b) high temperature 
(T = O.lhcJo/kB) and symmetrically biased voltage, (c) low 
temperature and asymmetric biased voltage {Vl — Vhias and 
Vr = 0). Other parameters include so = l.Ohwo, sd = lO/icJo 
for all plots. 


Besides the peak distances and peak width discussed 
above, other aspects characterizing the sidebands include 
the weights of the zero-phonon band and the number 
of peaks. The investigation of these properties mainly 
focuses on the effects of EPI strength, Fermi energy of 
the molecule, chemical potential and temperature of the 
leads. In general, the higher order peaks will be sup¬ 
pressed by the Frank-Condon factorIn the 
framework of NEGF, Ghen et al. found that the weight 
from zero-phonon band will decrease monotonically with 
increase of EPI strength and temperature while the 
weights of higher order sidebands will increase and then 
decrease^^. The chemical potentials of the leads will in¬ 
fluence the presence of the sidebands at both sides of 
the zero phonon peaks [Fig. 9, panel (a) and (b)]. If 
one keeps the chemical potential of one lead fixed and 
increases the chemical potential of the other lead, the 
phonon sidebands will appear only at one side of the 0- 
th order peak [Fig 9, panel (c)]. However, if one fixes 
the Fermi-level of the molecule and changes the chemical 
potentials of both leads, phonon sidebands will appear at 
both sides^^’^^. We note that, in this section, Fermi-level 
of the molecule is defined as (/i^ + /ii?)/2. 

In the framework of QME formalism described earlier, 
which is exact in the weak system-lead coupling limit, 
we find phonon sidebands for the single electronic level 
interacting with a single phonon mode. In this case the 
zero-phonon peak will occur at the renormalized resonant 
level due to polaron shift (^o ~ X^/huJo) and the side¬ 
bands will appear at every distance of hujQ. The peaks 
will appear at each side of zero-phonon peak under sym¬ 
metric change of the lead chemical potentials while only 
appear at one side if we fix chemical potential of one 
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lead. The EPI strength will not only shift the peaks, but 
also modulate the weights of the peaks. We find that 
upon increasing either EPI strength or temperature, the 
weight of the zero-phonon peak will decrease, which is 
consistent with the previous work^^. However, interest¬ 
ing phenomena happen when the renormalized level of 
the quantum dot gets close to the Eermi-level of the dot, 
[^0 — A^/cJo ~ (/^L + /^i?)/2], i.e., additional peaks ap¬ 
pear at each side. Those peaks have distance with 
the zero-order peak at the opposite side. However, the 
two major peaks really merge together, and those addi¬ 
tional peaks disappear again. In the case of asymmetric 
change of lead chemical potential (right most panel of 
Eig. 9 (c)), we found peaks appear at both sides when 
the zero-phonon peaks merge together. 

Another interesting perspective of the I — V charac¬ 
teristics is the phenomenon of the negative differential 
resistance (NDR), where the current decreases with the 
increase of voltage bias. The NEGE formalism predicts 
that NDR is impossible in ballistic electronic transport, 
but it will emerge in the presence of EPI. NDR has been 
both theoretically investigated^^’and experimen¬ 
tally measured^^^. An important reason of NDR is due 
to the redistribution of the molecular states. As discussed 
in the previous section, the zero-phonon band carries ma¬ 
jor portion of electronic current. The probability for the 
molecule to be in that state will be related to the chem¬ 
ical potential of the leads. If one increases the bias via 
increasing the chemical potential of one lead, one actually 
lifts the Eermi-level of the molecule as well. If one brings 
the Eermi-level of molecule far away from the eigenenergy 
of zero phonon state, the probability of the molecule to 
be in that state will decrease and thus the current will 
decrease. Based on this analysis one can draw several im¬ 
mediate conclusions: I. If one increases the bias simulta¬ 
neously for both leads in pace and keeps the Eermi-level 
of the molecule fixed, there will be no NDR. 2. If one 
treats the phonon fixed in the equilibrium distribution, 
NDR will not appear3. If the chemical-potential- 
varying lead couples stronger to the molecule than the 
chemical-potential-fixed lead, the redistribution will be 
more sensitive, thus the NDR will be enhanced. 

Eigure 10 shows the NDR predicted in the QME for¬ 
malism. We find that the NDR appears in the symmetric 
system-lead coupling. Moreover, NDR will be enhanced 
if the chemical-potential-varying lead couples stronger to 
the molecule than the other lead, but it will disappear in 
the other way around. We also find that NDR is most 
pronounced in the moderate EPI regime, while it is less 
significant in both weak and strong EPI regimes^^^. 


C. Vibrational effect on thermoelectricity 

In this part, we study the effects of EPI on the thermo¬ 
electric properties of the nano-conductor. We first look at 
the thermoelectric current, which is the electronic current 
induced by a temperature difference between the leads. 



VuasieV) VMasieV) VMasieV) 


FIG. 10. Prediction of NDR using QME formalism with the 
coupling strength (a) tjl = O.Iryi^, (b) r]L = Vr and (c) tjr = 
O.I?7l. For all plots Vr is fixed at 0 and Vl = Was- The 
temperature is fixed at T = 0.02huJo/kB for all leads. Other 
parameters are the same as Fig. 9 


Thermoelectric current exhibits quite different features 
comparing with voltage-bias current. It will increase 
monotonically and smoothly with the increasing of the 
temperature bias. Therefore, there will be no phonon 
sidebands. This is expected, because under temperature 
bias, the tunneling channel of the electron is always re¬ 
stricted to the molecular state that is close to the chemi¬ 
cal potential of the leads. The increasing of temperature 
will only excite more conducting electrons, but not be 
able to extend extra tunneling channels. Therefore, there 
will be no sudden change of thermoelectric current. Due 
to the same reason, there will be no NDR effect as the 
increasing of the temperature bias will always make more 
electrons to be involved in tunneling. The restriction on 
tunneling channels also makes the thermoelectric current 
much smaller than the voltage-bias current. 

The dependence of electronic conductance on the 
chemical potentials of the leads has been discussed in 
the Subsec. HIB. The sign of the Seebeck coefficient in¬ 
dicates that the currents will change direction with differ¬ 
ent chemical potential. Another interesting perspective 
is to find the dependence of currents on EPI strength. 
Eigure II shows the plot of voltage-bias current [panel 
(a)] and thermoelectric current [panel (b)] with the EPI 
strength A. Eor voltage-bias current, we find that the 
current decays with EPI strength in general. More pre¬ 
cisely, the maximum current one can achieve via adjust¬ 
ing the So is decreasing monotonically with EPI strength. 
This is due to the Erank-Gondon blockage of the current. 
However, for each we can find the enhancement of cur¬ 
rent due to EPI, and that is mainly because of the po- 
laron shift which can bring the electron resonance level 
into the conduction band from outside. Eor the ther¬ 
moelectric current, the profile is quite different. We see 
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FIG. 11. The dependence of voltage-bias current (a) and ther¬ 
moelectric current (b) on the EPI strength under different 
£o- For panel (a), jiL = hcjo, fJ^R = —hujo and = 

0.02hujo. For panel (b), — O.OShuJo, = 0.02hujo and 

fJ^L = fJ^R = 0 . 


that the current can change sign with the increase of EPI 
strength when sq is higher than the Fermi level, which 
indicates that the EPI can switch the charge carriers of 
the quantum dot between electrons and holes^^. For each 
eo there exist two optimized values of EPI strength such 
that the thermoelectric current maximizes. This opti¬ 
mized A shifts left with decrease of Gq until disappear 
one by one at the A = 0 end. 

One important quantity to describe the efficiency of 
the thermoelectric material is the figure of merits Z^T. 
It is related to the electronic conductance Gg, Seebeck 
coefficient S', thermal conductance Gh via the formula 
ZeT = GeS‘^T/Gh‘ Here we use the notation ZqT 
to denote the figure of merits of the system by ignor¬ 
ing the thermal conductance due to phonons. So Gh 
here only takes account the thermal conductance due 
to electrons. The effect of phonons on the figure of 
merit ZqT is rather complicated, which is closely related 
to the electron Fermi energy^^’^^^’^^^’^^^, the phonon 
energythe temperature and chemical potentials of 
the leads^^^’^^^’^^^’^^^. Figure 12 shows the dependence 
of the electronic conductance, thermal conductance, See¬ 
beck coefficient and figure of merit on the electron en¬ 
ergy So and EPI strength A. The major effect of EPI is 
on the thermal conductance of the molecular junction^^^. 
The EPI can open extra channels, from which high en¬ 
ergy electrons can tunnel from hotter lead to colder lead, 
while low energy electrons tunnel in the opposite direc¬ 
tion. Therefore, the thermal conductance is enhanced 
while the electronic conductance is not affected too much. 
Due to the increase of the thermal conductance, Z^T will 
be suppressed drastically^^^’^^^. Though the figure of 


FIG. 12. The contour plot of electronic conductance Ge, ther¬ 
mal conductance Gh, Seebeck coefficient S and figure of merit 
ZeT with respect to Sq and A. The parameters are: /x = 0, 
T = bhuJo/kB- For this plot, the phonon mode is coupled 
to its environments at temperature T with coupling strength 
TjE = 0.1?7. Figure adapted with permission from Phys. Rev. 
B, 91, 045410 (2015). Gopyrighted by the American Physical 
Society. 


merits Z^T will be reduced quickly under the influence 
of phonon scattering in weak EPI regime, it will gradu¬ 
ally saturate at strong EPI regime. We also would like 
to point out that when the electron energy is close to the 
Fermi energy of the quantum dot, the Seebeck coefficient 
and hence the Z^T will become very small. However, 
the phonon scattering can renormalized the electron en¬ 
ergy via polaron shift and thus the Seebeck coefficient 
can be enhanced. As a result, EPI can enhance Z^T in 
this particular parameter regime. 


D. Local heating 

Local heating is an important phenomenon in molec¬ 
ular junction, not only due to its own importance af¬ 
fecting the stability of the system, but also due to its 
close relation to the phonon sidebands^^^, NDR and ther¬ 
moelectric effect^^’^^^. The distribution of the phonon 
states in current-carrying system can be far away from 
equilibrium^^^, or in some cases, may even lead to phonon 
instability^^’'^^’^^^’^^^. Phonons can be excited signif¬ 
icantly by the voltage bias^^^, and in turn affect the 
I — V characteristics. At each peak of the phonon side¬ 
bands, one can actually find a vibrational excitation 
eventPrevious study also found that the local 
heating can enhance the thermoelectric efficiency^^’^^^. 
However, in most cases local heating is not preferred 
because it will affect the stability of the system^"^ and 
introduce noise to the measurements^^^”^^^. Therefore, 
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lot of effort has been put into cooling the system using 
electronic current, such as using superconducting single 
electron transistor^^^, or double quantum dots^^^. 

Here we study the effects of the electronic current on 
the phonon mode of the nano-conductor. To investigate 
the local heating, effective temperatures are usually de¬ 
fined in various In this part, since we 

only have one single phonon mode, instead of defining the 
effective temperature, we use average phonon number to 
characterize the local heating effect. The way to spec¬ 
ify the electronic current induced heating is to compare 
the nonequilibrium phonon number rineq with equilib¬ 
rium phonon number When the molecule is weakly 
coupled to the leads, the molecular states statistics will 
follow canonical distribution p = e“^^‘^/Tr(e“^^‘^) and 
thus the equilibrium phonon number can be calculated 
exactly as^^ 


1 


^ea — 


+ 


\^/{hujo) 


^PhujQ _ I ^I3 {£o-X^/( huJo)) ^ l' 


(75) 


The first term is the Bose-Einstein distribution function, 
the second term is a correction due to the polaron en¬ 
ergy shift. The nonequilibrium phonon number can be 
calculated from the nonequilibrium reduced density ma¬ 
trix obtained from the QME. Eigure 13 shows the dif¬ 
ference of phonon numbers under voltage bias (top) and 
temperature bias (down). Eor voltage bias. An is always 
positive which indicates that the system is always heated 
up. The yellow regime where An ~ 0 is the regime where 
the electronic current vanishes. When there is electronic 
current passing through, the heating effect is generally 
more pronounced in stronger EPI regime. However, for 
the temperature bias, we find both heating (An > 0) and 
cooling regimes (An < 0). Therefore, the local heating 
effect is not only related to the magnitude of electronic 
current, but also related to the energy each electron car¬ 
ries when it tunnels into the molecule. Eor the thermo¬ 
electric current, low energy electron can tunnel from the 
cooler lead to the molecule, absorb a phonon and tunnel 
to the hotter lead, resulting in cooling of the molecule. 
Such process is impossible in the voltage-bias case in the 
present setup as the electron is always flowing from the 
higher chemical potential side to the lower chemical po¬ 
tential side. 
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FIG. 13. The correction of phonon number due to elec¬ 
tronic current under voltage bias (a) and temperature bias 
(b). The parameters are: /xl = —pn = 2huJo, = T + AT, 
T^ = T — AT, AT = Shojo/ks and T = bhujo/kB ■ The 
phonon is weakly coupled to its own environment at temper¬ 
ature T. Figure adapted with permission from Phys. Rev. 
B, 91, 045410 (2015). Copyrighted by the American Physical 
Society. 


leads. The Langevin equation applies to the weak EPI 
regime, since similar to Sec. HI our expansion parame¬ 
ter is the interaction matrix M. However, the derivation 
is not limited to the form of i^tot- Hollowing the same 
procedure we discuss below, we can also do an adiabatic 
expansion of the electron influence functional over the 
velocity of the ions. In that case, the Langevin equation 
applies to slow ions, but the EPI could be of arbitrary 
magnitude^^“^^’^^^ . 


V. CURRENT-INDUCED SEMI-CLASSICAL 
LANGEVIN DYNAMICS 

In the previous two sections, we mainly look at the ef¬ 
fect of phonons on the electric and thermoelectric trans¬ 
port properties of electrons, which is also the focus of 
most published work. But, to study the current-induced 
forces, and their effect on atomic dynamics, we need to 
turn around. In this section, starting from the total 
Hamiltonian iLtot, we derive a semi-classical Langevin 
equation to describe the atomic dynamics of the system, 
coupled to both phonon and (nonequilibrium) electron 


The advantage of the Langevin equation approach is 
that, we can easily include the anharmonic phonon in¬ 
teraction, as in other molecular dynamics method. The 
anharmonic interaction is crucial in dealing with current- 
induced dynamics. This is because the phonon modes 
that interact with electrons are normally high-frequency 
ones, while the low frequency modes conduct heat from 
the system to surrounding electrodes more effectively. 
The energy transfer from high to low frequency modes is 
possible only when anharmonic interaction is included. 
Although possible, it is not a trivial task to incorporate 
anharmonic interactions in the NEGE or QME approach. 
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A. Initial states and reduced density matrices 


We assume at a remote past to and earlier time, the 
central system is decoupled with the leads and there is 
also no EPI, so that the electrons and phonons are also 
decoupled. The density matrix is assumed to be product 
of known equilibrium states. For example, the left lead 
is specified by 




(76) 

(77) 


Exactly what to take for the center is not important as 
for steady state with to — oo, the results do not depend 
on it (except maybe in very subtle cases). 

The density matrix (of the whole system) is governed 
by the von Neumann equation, and formally we can write 


where 


p(t) = U{t,to)p{to)U{to,t), 


U{t, to) = -^‘o 


(78) 


(79) 


assuming t > to- For the other case of t < to, the time- 
order operator T should be replaced by the anti-time 
order operator. We are interested only in the center, 
so the leads degrees of freedom will be traced out. For 
notational simplicity, we assume only one left lead. The 
result for two or more leads is trivially generalized. We 
define 


p{t) = TTLp{t), (80) 

Pp(0 = Trep(0, (81) 

where the first reduced density matrix only eliminates the 
lead, while the last one eliminates the electrons as well, 
leaving only an effective density matrix for the phonons. 
The procedure to eliminate the lead for both the electrons 
and phonons follows the standard method of Feynman 
and Vernon^^^, except that we need to be careful for the 
electrons which are fermions. Since there is no coupling 
between electrons and phonons in the leads, the phonon 
and electron degrees of freedom can be done separately 
(the initial density matrix is a product of the two). 


B. Influence functional for phonons 

The matrix elements of the density operator p(t) is 
taken in the basis of the coordinates u. Following the 
standard treatmentthe density matrix is then 
given by 

{u'\p{t)\u) (X. j pq{u'o,uo). (82) 

The path C 2 consists of two segments (see Fig. 14(a)) 
running from time t with coordinate variable u back to 
to at variable uo^ and then second segment running from 


(a) a 

I I I I l >l I I H-0 _^ 

I I I I I <1 I I I O 

(b) c 

I I I I l>l I _ 

I I I I M l I 11 ^ ^ 

to t 

FIG. 14. Two types of paths in the Feynman-path integrals, 
(a) two-segment path C 2 . (b) Keldysh type closed contour C. 
The ticks represent the discretized integration variables; the 
open circles are not integrated. 


to with variable u'o forward to time t with variable u'. 
Since the arbitrary constant in the proportionality can 
be fixed by normalization (Trp = 1), we need not specify 
precisely the measure associated with V[u]. The path 
integral integrates all the intermediate variables except 
the two open ends at time t. 

To eliminate the lead, the phonon Lagrangian is split 
into jC = jC^ -i- jC^ — V where V is the coupling poten¬ 
tial energy term between the lead and center, given by 
V = . The integration volume elements are 

also split V[u] = V[u^]V[u^]. The initial distribution 
is assumed product type Po ^ Po • Taking the trace of 
Eq. (82) means that we identify and {u')^ as the same 
variable and integrate it out. As far as variable 
is concerned, the path is of Keldysh type, i.e., running 
from to to t from above and then back from t to to, as 
shown in Fig. 14(b). The reduced density matrix is then 


~p{t)= f 

(83) 

with the influence functional 


(84) 

The above expression can be further simplified. Firstly, 
the initial distribution of the lead is in thermal equilib¬ 
rium, po oc . Secondly, we can rewrite the path 

integral formula back into the operator form; thus, we 
obtain 


I[u^{t)\ = Tr 


e „ -{i/h)J^{H^+v)dT 


Zl 


-Tee 




(85) 


where Zl the canonical partition function, Tc is con¬ 
tour order operator, the subscript eq.L stands for equi¬ 
librium average with respect to the left lead. One more 
transformation can be made to simplify it further. The 
time-dependence (or rather, contour time r dependence) 
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is understood to be in the Heisenberg picture governed 
by the Hamiltonian + V (r) which is different in the 
forward and backward direction. We can work in the 
interaction picture, thus eliminate the explicit from 
the formula. The resulting equation is 

I[u^it)] = fo (86) 

This is the same equation [Eq. (5)] given in Ref. 148. 


The inffuence functional can be calculated explicitly 
since the interaction H is a quadratic form, and the con¬ 
tour operator naturally leads to contour ordered Green’s 
functions and Wick’s theorem is valid. We define the 
contour ordered phonon Green’s function of the lead as 

dL{T,T') = -^{TcU^{T)u^{T')'^)eq.h- (87) 

Expanding the exponential (or a cumulant expansion, ex¬ 
panding and taking logarithm), we get 


I[u'^{t)] = {Tc{ 1- ^ f ^i(T)dT-^ f f dT'Vl{T)Vl{T') + ■■■)) 
J C <7 C J c 

= Tc(l - T ^ ^ M^(r)^y^^(Tc«^(r)M^(rO^)V^^M^(r')drdT' + ■ ■ ■) 

= Tc(l-“ j j dT'u'^{T)'^IlL{T,T')u'^{T') -) 


The first order term (and all odd order in V terms) van¬ 
ishes, because {u^) = 0. We have defined the ubiquitous 
contour ordered self-energy due to the lead as 

ni(T, r') = V^^dLir, (89) 

It can be shown that the last line in the above derivation 
is an exact result. 

It is instructive to clarify and compare with other no¬ 
tations used in the literature. A commonly used notation 
is (e.g.. Ref. 103) 

1 T 

In/=-- / dti / dt 2 (h) - u~(ti)) 

d Jto Jto 

X [L{ti,t2)u^{t2) - L^{ti,t2)u~{t2)] . (90) 

Using the rules dr crdt^ and n(r, r') 

H^’^ u^{r) u^{t) where a = -f or — for the 

upper and lower branch, the relations among the var¬ 
ious self-energies, and symmetry relation n^(t,t') = 
n^^(t', t), we can rewrite Eq. (88) in the form of Eq. (90). 
By comparison, we find 

L(t, d) = m>(t, t'), T*(t, d) = iii< (t, d). (91) 

a{t,d) = L(t^d) is the notation used by Schmid^^^. 

C. Influence functional for electrons 

The derivation of the influence functional for the 
electrons is similar except that we have to deal with 
grassmann integrals^^^”^^^. We follow the approach of 
Weinberg^^^. To do the trace over the leads we need a 


specific representations for the operator p. Eor the elec¬ 
trons, we use the coherent state characterized by grass¬ 
mann numbers such that 


T7 

II 

(92) 

The hat denotes operator, without hat, it is a grassmann 
number. The state c) has explicit form given as 


(93) 

i 

(94) 

The orthogonality is in the form 


(c'ic)= n(ci -c'). 

i 

(95) 

Similar states for the creation operator ct can 

be con- 

structed with eigenvalue (a grassmann number) c, given 
the following results (similar to inner product of eigen 
states of u and its conjugate momentum p and complete¬ 
ness of the eigenstates). 

(cc) = ---e 

(96) 

{c\c) = , 

(97) 

1 = / \c)tli-dcj){c\, 

(98) 

1 = / \c)'[lidcj){c\. 

(99) 


The • • • is an extra + or — sign factor which we’ll not 
keep track. The tilde over the product sign means that 
order is in exactly the opposite canonical order [e.g., that 
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of Eq. (95)]. With the above very sketching outline, the 
fermion evolution operator 

U{t, to) = -^‘0 (100) 

can be represented as a path integral of the form 

j (101) 

with the action 

Se = — J dr He — . (102) 

The lead influence functional can be then obtained with 
the same procedure as that for the phonons. The result 
involves an integral kernel which is exactly the contour 
ordered self-energy of the lead ^^(r, r'). 


D. Reduced density matrix for phonon in center 

Putting things together, the reduced density matrix, 
when the leads are eliminated, has the form 

J V[u^, (p, c^Y^I^po{{PoT■, reF■, (co)*^; Co*^, Y) 

where 

S = Sp Sqp + 5'p + S'e, (103) 

= j dr u - Vn{u)^ , (104) 

= J dr (^hc^ ^ , (105) 

5'^ = - f dr'^^UkC^M^c, (106) 

dC2 ^ 

= [ dr f dr'u{r)^ULir,H)u(r'), (107) 

si = - ! dr f dr'c(r)^EL(r,r')c(r'). (108) 

JC2 Jc2 

The ordinary number (column vector) u and grassmann 
number c, c involve only the degrees of freedom of the 
center. For notational simplicity, we have dropped the 
superscript C. Note that the electron terms do not have 
the characteristic factor 1/2 as c and c are independent 
variables. 

The dependence on c and c is a bi-linear form, thus the 
path integral over them can be done analytically. This 
gives the reduced density matrix of the phonon only, as 

{u'WYW) ^ j (109) 


The influence functional to the phonons due to electrons 
is given by 

Ip oc det (^{r,r') ^lih^ - - Uk{r)M^^ - T^L{r,r') 

( 110 ) 

Interpreting the r in the above as Keldysh variable de¬ 
fined on C has a problem. As agreed, the contour is 
supposed to be on C 2 with the to end connected with the 
initial distribution of the electrons po at the center. How¬ 
ever, if we assume that in the limit to ^ — oc the results 
should not depend on the distribution of the center, we 
can ignore this initial distribution and it is completely 
fixed by the lead. But we cannot give a mathematically 
sound justification here. 

Similar to that in Ref. 154 for the field theory of quan¬ 
tum electrodynamics, we want to put the influence func¬ 
tional in an exponential form. This can be done using 
the formula, Det (A) = expansion of the 

function ln(l + x) = x — x^/2 + • • • for small x, given 


n+l 


Ip = det(Go ^+ 2 /) = det(Go ^) exp (-Tr[(Go2/)"] 


Vn=l 


( 111 ) 


where 


Go = S{t,t') 


-5]L(r,r'), (112) 


y = -S{t, t') Ufe(r)M''. 


(113) 


Gq ^ and y are matrices indexed by lattice sites j as well 
as contour time r. And, if the proper metric for a dis¬ 
cretization of the time is chosen so that det(G^^) can be 
meaningful, we can identify Go as the electron contour 
ordered Green’s function when there is no EPI defined 
in Sec. III. Since det(Go ^) is independent of the ef¬ 
fective action for the phonon is only determined by the 
exponential factor, which is a polynomial (functional) in 
u. With some caveat regarding the initial distribution, 
Eqs. (109) and (111) offer a formally exact solution to 
the problem. 

The first two terms, written out explicitly in terms of 
the contour ordered Green’s functions are 


T^lGoy] = [Go{t,t)M’^] Uk{T), (114) 

k 

-^Tr [(Go?/)^] =-^ E JdTdT'dT"dT"'Goji{T,T') 

VlrnY .r'')Go^YT\T''')yrni{T"', t) 

= -^ / J dT'u{T)'^Ilep{T,T')u{T'), (115) 


with 


nep'(Gr') = -ihTr \Go{t',t)M'^Go{t,t')M'^'] . (116) 
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E. Semi-classical approximation 


If we ignore the linear term in u which produces a con¬ 
stant force, the effect of which is to shift the equilibrium 
positions, and also neglect higher order contributions, we 
end up with a quadratic form for the effective action 

S'eff = J dr - ^u'^K^u - Vn{u)^ 

J C2 C2 

with Iltot = 11^ + + Hep as in Eq. (20). We have also 

included the right phonon lead. A generalized Langevin 
equation can be derived from the above action^^d49 


il = —K^u + Fn 





(118) 


the electrons^^. Moreover, the power of the Langevin 
approach is to be able to include the anharmonic 
phonon-phonon interactions classically, and treat the 
EPI quantum-mechanically. This enables one to study 
the energy transport between different phonon modes, 
between electrons and phonons at the same time. The 
exploration of its power and potential is still under way. 

As an example, we consider the heat generation in 
a 4 X 2 graphene armchair ribbon, see Ref. 157 for 
the definition of structure parameters. The electronic 
structure and EPI matrix are obtained from a combined 
SIESTA^^^ -h TranSIESTA^^^ -1- Inelastica^^ calculation, 
while the Brenner potential is used for the inter-atomic 
interaction. To reduce the simulation time, we have ig¬ 
nored the energy-dependence of the electronic structure. 
As a result, the electronic friction becomes time-local. 
The Langevin equation becomes (after an integration by 
part) 


where Fn = —dVnjdu^ is the retarded total self¬ 
energy, and the noises satisfy 

m) = 0 , ( 119 ) 

= id- -1')+ n<t(t - 1')) 

= (120) 

We note that the effect of the electron leads to the 
phonons has exactly the same form as that of the phonon 
leads. The self-energy consists of a sum of contributions 
of the two sources. 


F. Applications 

Before discussing the applications, we note that 
similar generalized Langevin equation as Eq. (118) 
can be derived by doing an adiabatic expan¬ 
sion over the momenta of the ions to the 2nd 
order^^~^^’^^^. These equations have been used in 
different perspective^^“^^’^^’^^^’^^^’^^^~^^^. Its most 
important feature is the inclusion of the quantum nature 
of the electron and phonon leads. Eor example, the 
zero point motions of atoms are correctly taken into 
account, and proved critical in determining the thermal 
and structural properties of materials made from light 
elements^^^”^^^’^^^’^^^. Including the correct Bose 
distribution of the phonons opens a way to study the 
quantum ballistic phonon transport by doing classical 
molecular dynamics. In Refs. 156 and 157 the transition 
from ballistic to diffusive phonon thermal transport 
is studied using this approach. In Refs. 17 and 34, 
including the nonequilibrium electrons, current-induced 
dynamics have been studied. Several interesting effects 
have been predicted or confirmed, and their effects on 
the stability of the system are studied. Eor example, 
it has been shown that (1) the current-induced forces 
are not conservative^^, (2) the atoms feel an effective 
magnetic force, originating from the Berry phase of 


u = F^ 



r{t—t')u{t')dt' — hr]u—eV^ (121) 


where F^ is the force from the second-generation Bren¬ 
ner potential, dT{t)/dt = 11^ (t) is the phonon retarded 
self-energy due to two leads, V is applied bias voltage. 
The eV^~u term gives a nonconservative force as is 
antisymmetric. ^ F is the noise due to left and 
right phonon leads, while / is the noise due to electron 
bath. The expression for the electronic friction and noise 
correlation is the same as Eqs. (56-63) in Ref. 18. Eur- 
ther implementation details can be found in Ref. 157. 
The phonon heat current is calculated using 




t')u{t')dt' F 


( 122 ) 


where a = L^R. In steady state, the energy flow bal¬ 
ances, and the heat generation is calculated according to 
Q = - Jp«. 

The rate of heat generation is plotted in Eig. 15. The 
result is for the same configuration as shown in Ref. 157 
of Eig. 4. Each data point takes about 4 days on a typ¬ 
ical Opteron CPU. The error bars are quite small. The 
heat generation at zero bias should be zero. However, we 
get a small value. This has to do with the cut-off used 
in the noise for the electrons. We have used an abrupt 
cut-off for the spectrum at huj = 1.29 eV. The calculation 
demonstrates the feasibility of computing the Joule heat¬ 
ing current, intrinsically a quantum effect at nanoscale, 
by classical molecular dynamics. 


VI. CONCLUSIONS 

In summary, using a tight-binding-like Hamiltonian for 
the EPI, Eq. (8), we have introduced three different ap¬ 
proaches to study the effect of EPI in different parame¬ 
ter regimes. We focused on the electronic, phononic, and 
thermoelectric transport properties of nano-conductors. 
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lems are: (1) application of the NEGF and QME ap¬ 
proach to nonequilibrium thermoelectric transport to 
study the thermoelectric efficiency at finite power out¬ 
put, (2) application of the QME approach to look at 
system where both EPI and electron-electron interaction 
are important, (3) combining the generalized Langevin 
equation with first-principles or tight-binding electronic 
structure package to study current-induced dynamics in 
realistic nano-conductors, especially to explore how the 
electron-dissipated heat is transferred in and out of the 
nano-conductor. With these available tools, more inter¬ 
esting and important systems can be investigated. 
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